<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>slaqtr.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="SLAQTR.1"></a><a href="slaqtr.f.html#SLAQTR.1">SLAQTR</a>( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
     $                   INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK auxiliary routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      LOGICAL            LREAL, LTRAN
      INTEGER            INFO, LDT, N
      REAL               SCALE, W
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      REAL               B( * ), T( LDT, * ), WORK( * ), X( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="SLAQTR.20"></a><a href="slaqtr.f.html#SLAQTR.1">SLAQTR</a> solves the real quasi-triangular system
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">               op(T)*p = scale*c,               if LREAL = .TRUE.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  or the complex quasi-triangular systems
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">             op(T + iB)*(p+iq) = scale*(c+id),  if LREAL = .FALSE.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  in real arithmetic, where T is upper quasi-triangular.
</span><span class="comment">*</span><span class="comment">  If LREAL = .FALSE., then the first diagonal block of T must be
</span><span class="comment">*</span><span class="comment">  1 by 1, B is the specially structured matrix
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 B = [ b(1) b(2) ... b(n) ]
</span><span class="comment">*</span><span class="comment">                     [       w            ]
</span><span class="comment">*</span><span class="comment">                     [           w        ]
</span><span class="comment">*</span><span class="comment">                     [              .     ]
</span><span class="comment">*</span><span class="comment">                     [                 w  ]
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  op(A) = A or A', A' denotes the conjugate transpose of
</span><span class="comment">*</span><span class="comment">  matrix A.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  On input, X = [ c ].  On output, X = [ p ].
</span><span class="comment">*</span><span class="comment">                [ d ]                  [ q ]
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  This subroutine is designed for the condition number estimation
</span><span class="comment">*</span><span class="comment">  in routine <a name="STRSNA.45"></a><a href="strsna.f.html#STRSNA.1">STRSNA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LTRAN   (input) LOGICAL
</span><span class="comment">*</span><span class="comment">          On entry, LTRAN specifies the option of conjugate transpose:
</span><span class="comment">*</span><span class="comment">             = .FALSE.,    op(T+i*B) = T+i*B,
</span><span class="comment">*</span><span class="comment">             = .TRUE.,     op(T+i*B) = (T+i*B)'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LREAL   (input) LOGICAL
</span><span class="comment">*</span><span class="comment">          On entry, LREAL specifies the input matrix structure:
</span><span class="comment">*</span><span class="comment">             = .FALSE.,    the input is complex
</span><span class="comment">*</span><span class="comment">             = .TRUE.,     the input is real
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          On entry, N specifies the order of T+i*B. N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  T       (input) REAL array, dimension (LDT,N)
</span><span class="comment">*</span><span class="comment">          On entry, T contains a matrix in Schur canonical form.
</span><span class="comment">*</span><span class="comment">          If LREAL = .FALSE., then the first diagonal block of T must
</span><span class="comment">*</span><span class="comment">          be 1 by 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDT     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the matrix T. LDT &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  B       (input) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, B contains the elements to form the matrix
</span><span class="comment">*</span><span class="comment">          B as described above.
</span><span class="comment">*</span><span class="comment">          If LREAL = .TRUE., B is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  W       (input) REAL
</span><span class="comment">*</span><span class="comment">          On entry, W is the diagonal element of the matrix B.
</span><span class="comment">*</span><span class="comment">          If LREAL = .TRUE., W is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  SCALE   (output) REAL
</span><span class="comment">*</span><span class="comment">          On exit, SCALE is the scale factor.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  X       (input/output) REAL array, dimension (2*N)
</span><span class="comment">*</span><span class="comment">          On entry, X contains the right hand side of the system.
</span><span class="comment">*</span><span class="comment">          On exit, X is overwritten by the solution.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          On exit, INFO is set to
</span><span class="comment">*</span><span class="comment">             0: successful exit.
</span><span class="comment">*</span><span class="comment">               1: the some diagonal 1 by 1 block has been perturbed by
</span><span class="comment">*</span><span class="comment">                  a small number SMIN to keep nonsingularity.
</span><span class="comment">*</span><span class="comment">               2: the some diagonal 2 by 2 block has been perturbed by
</span><span class="comment">*</span><span class="comment">                  a small number in <a name="SLALN2.95"></a><a href="slaln2.f.html#SLALN2.1">SLALN2</a> to keep nonsingularity.
</span><span class="comment">*</span><span class="comment">          NOTE: In the interests of speed, this routine does not
</span><span class="comment">*</span><span class="comment">                check the inputs for errors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            NOTRAN
      INTEGER            I, IERR, J, J1, J2, JNEXT, K, N1, N2
      REAL               BIGNUM, EPS, REC, SCALOC, SI, SMIN, SMINW,
     $                   SMLNUM, SR, TJJ, TMP, XJ, XMAX, XNORM, Z
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      REAL               D( 2, 2 ), V( 2, 2 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            ISAMAX
      REAL               SASUM, SDOT, <a name="SLAMCH.116"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANGE.116"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>
      EXTERNAL           ISAMAX, SASUM, SDOT, <a name="SLAMCH.117"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANGE.117"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SAXPY, <a name="SLADIV.120"></a><a href="sladiv.f.html#SLADIV.1">SLADIV</a>, <a name="SLALN2.120"></a><a href="slaln2.f.html#SLALN2.1">SLALN2</a>, SSCAL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Do not test the input parameters for errors
</span><span class="comment">*</span><span class="comment">
</span>      NOTRAN = .NOT.LTRAN
      INFO = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set constants to control overflow
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.139"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SMLNUM = <a name="SLAMCH.140"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> ) / EPS
      BIGNUM = ONE / SMLNUM
<span class="comment">*</span><span class="comment">
</span>      XNORM = <a name="SLANGE.143"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>( <span class="string">'M'</span>, N, N, T, LDT, D )
      IF( .NOT.LREAL )
     $   XNORM = MAX( XNORM, ABS( W ), <a name="SLANGE.145"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>( <span class="string">'M'</span>, N, 1, B, N, D ) )
      SMIN = MAX( SMLNUM, EPS*XNORM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute 1-norm of each column of strictly upper triangular
</span><span class="comment">*</span><span class="comment">     part of T to control overflow in triangular solver.
</span><span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = ZERO
      DO 10 J = 2, N
         WORK( J ) = SASUM( J-1, T( 1, J ), 1 )
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      IF( .NOT.LREAL ) THEN
         DO 20 I = 2, N
            WORK( I ) = WORK( I ) + ABS( B( I ) )
   20    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      N2 = 2*N
      N1 = N
      IF( .NOT.LREAL )
     $   N1 = N2
      K = ISAMAX( N1, X, 1 )
      XMAX = ABS( X( K ) )
      SCALE = ONE
<span class="comment">*</span><span class="comment">
</span>      IF( XMAX.GT.BIGNUM ) THEN
         SCALE = BIGNUM / XMAX
         CALL SSCAL( N1, SCALE, X, 1 )
         XMAX = BIGNUM
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( LREAL ) THEN
<span class="comment">*</span><span class="comment">
</span>         IF( NOTRAN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve T*p = scale*c
</span><span class="comment">*</span><span class="comment">
</span>            JNEXT = N
            DO 30 J = N, 1, -1
               IF( J.GT.JNEXT )
     $            GO TO 30
               J1 = J
               J2 = J
               JNEXT = J - 1
               IF( J.GT.1 ) THEN
                  IF( T( J, J-1 ).NE.ZERO ) THEN
                     J1 = J - 1
                     JNEXT = J - 2
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( J1.EQ.J2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Meet 1 by 1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale to avoid overflow when computing
</span><span class="comment">*</span><span class="comment">                     x(j) = b(j)/T(j,j)
</span><span class="comment">*</span><span class="comment">
</span>                  XJ = ABS( X( J1 ) )
                  TJJ = ABS( T( J1, J1 ) )
                  TMP = T( J1, J1 )
                  IF( TJJ.LT.SMIN ) THEN
                     TMP = SMIN
                     TJJ = SMIN
                     INFO = 1
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( XJ.EQ.ZERO )
     $               GO TO 30
<span class="comment">*</span><span class="comment">
</span>                  IF( TJJ.LT.ONE ) THEN
                     IF( XJ.GT.BIGNUM*TJJ ) THEN
                        REC = ONE / XJ
                        CALL SSCAL( N, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
                  X( J1 ) = X( J1 ) / TMP
                  XJ = ABS( X( J1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale x if necessary to avoid overflow when adding a
</span><span class="comment">*</span><span class="comment">                 multiple of column j1 of T.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( XJ.GT.ONE ) THEN
                     REC = ONE / XJ
                     IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
                        CALL SSCAL( N, REC, X, 1 )
                        SCALE = SCALE*REC
                     END IF
                  END IF
                  IF( J1.GT.1 ) THEN
                     CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
                     K = ISAMAX( J1-1, X, 1 )
                     XMAX = ABS( X( K ) )
                  END IF
<span class="comment">*</span><span class="comment">
</span>               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Meet 2 by 2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Call 2 by 2 linear system solve, to take
</span><span class="comment">*</span><span class="comment">                 care of possible overflow by scaling factor.
</span><span class="comment">*</span><span class="comment">
</span>                  D( 1, 1 ) = X( J1 )
                  D( 2, 1 ) = X( J2 )
                  CALL <a name="SLALN2.251"></a><a href="slaln2.f.html#SLALN2.1">SLALN2</a>( .FALSE., 2, 1, SMIN, ONE, T( J1, J1 ),
     $                         LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
     $                         SCALOC, XNORM, IERR )
                  IF( IERR.NE.0 )
     $               INFO = 2
<span class="comment">*</span><span class="comment">
</span>                  IF( SCALOC.NE.ONE ) THEN
                     CALL SSCAL( N, SCALOC, X, 1 )
                     SCALE = SCALE*SCALOC
                  END IF
                  X( J1 ) = V( 1, 1 )
                  X( J2 ) = V( 2, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2))
</span><span class="comment">*</span><span class="comment">                 to avoid overflow in updating right-hand side.
</span><span class="comment">*</span><span class="comment">
</span>                  XJ = MAX( ABS( V( 1, 1 ) ), ABS( V( 2, 1 ) ) )
                  IF( XJ.GT.ONE ) THEN
                     REC = ONE / XJ
                     IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
     $                   ( BIGNUM-XMAX )*REC ) THEN
                        CALL SSCAL( N, REC, X, 1 )
                        SCALE = SCALE*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Update right-hand side
</span><span class="comment">*</span><span class="comment">
</span>                  IF( J1.GT.1 ) THEN
                     CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
                     CALL SAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
                     K = ISAMAX( J1-1, X, 1 )
                     XMAX = ABS( X( K ) )
                  END IF
<span class="comment">*</span><span class="comment">
</span>               END IF
<span class="comment">*</span><span class="comment">
</span>   30       CONTINUE
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve T'*p = scale*c
</span><span class="comment">*</span><span class="comment">
</span>            JNEXT = 1
            DO 40 J = 1, N
               IF( J.LT.JNEXT )
     $            GO TO 40
               J1 = J
               J2 = J
               JNEXT = J + 1
               IF( J.LT.N ) THEN
                  IF( T( J+1, J ).NE.ZERO ) THEN
                     J2 = J + 1
                     JNEXT = J + 2
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( J1.EQ.J2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 1 by 1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale if necessary to avoid overflow in forming the
</span><span class="comment">*</span><span class="comment">                 right-hand side element by inner product.
</span><span class="comment">*</span><span class="comment">
</span>                  XJ = ABS( X( J1 ) )
                  IF( XMAX.GT.ONE ) THEN
                     REC = ONE / XMAX
                     IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
                        CALL SSCAL( N, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  X( J1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X, 1 )
<span class="comment">*</span><span class="comment">
</span>                  XJ = ABS( X( J1 ) )
                  TJJ = ABS( T( J1, J1 ) )
                  TMP = T( J1, J1 )
                  IF( TJJ.LT.SMIN ) THEN
                     TMP = SMIN
                     TJJ = SMIN
                     INFO = 1
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( TJJ.LT.ONE ) THEN
                     IF( XJ.GT.BIGNUM*TJJ ) THEN
                        REC = ONE / XJ
                        CALL SSCAL( N, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
                  X( J1 ) = X( J1 ) / TMP
                  XMAX = MAX( XMAX, ABS( X( J1 ) ) )
<span class="comment">*</span><span class="comment">
</span>               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 2 by 2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale if necessary to avoid overflow in forming the
</span><span class="comment">*</span><span class="comment">                 right-hand side elements by inner product.
</span><span class="comment">*</span><span class="comment">
</span>                  XJ = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ) )
                  IF( XMAX.GT.ONE ) THEN
                     REC = ONE / XMAX
                     IF( MAX( WORK( J2 ), WORK( J1 ) ).GT.( BIGNUM-XJ )*
     $                   REC ) THEN
                        CALL SSCAL( N, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  D( 1, 1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X,
     $                        1 )
                  D( 2, 1 ) = X( J2 ) - SDOT( J1-1, T( 1, J2 ), 1, X,
     $                        1 )
<span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SLALN2.370"></a><a href="slaln2.f.html#SLALN2.1">SLALN2</a>( .TRUE., 2, 1, SMIN, ONE, T( J1, J1 ),
     $                         LDT, ONE, ONE, D, 2, ZERO, ZERO, V, 2,
     $                         SCALOC, XNORM, IERR )
                  IF( IERR.NE.0 )
     $               INFO = 2
<span class="comment">*</span><span class="comment">
</span>                  IF( SCALOC.NE.ONE ) THEN
                     CALL SSCAL( N, SCALOC, X, 1 )
                     SCALE = SCALE*SCALOC
                  END IF
                  X( J1 ) = V( 1, 1 )
                  X( J2 ) = V( 2, 1 )
                  XMAX = MAX( ABS( X( J1 ) ), ABS( X( J2 ) ), XMAX )
<span class="comment">*</span><span class="comment">
</span>               END IF
   40       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span>         SMINW = MAX( EPS*ABS( W ), SMIN )
         IF( NOTRAN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve (T + iB)*(p+iq) = c+id
</span><span class="comment">*</span><span class="comment">
</span>            JNEXT = N
            DO 70 J = N, 1, -1
               IF( J.GT.JNEXT )
     $            GO TO 70
               J1 = J
               J2 = J
               JNEXT = J - 1
               IF( J.GT.1 ) THEN
                  IF( T( J, J-1 ).NE.ZERO ) THEN
                     J1 = J - 1
                     JNEXT = J - 2
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( J1.EQ.J2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 1 by 1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale if necessary to avoid overflow in division
</span><span class="comment">*</span><span class="comment">
</span>                  Z = W
                  IF( J1.EQ.1 )
     $               Z = B( 1 )
                  XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
                  TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
                  TMP = T( J1, J1 )
                  IF( TJJ.LT.SMINW ) THEN
                     TMP = SMINW
                     TJJ = SMINW
                     INFO = 1
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( XJ.EQ.ZERO )
     $               GO TO 70
<span class="comment">*</span><span class="comment">
</span>                  IF( TJJ.LT.ONE ) THEN
                     IF( XJ.GT.BIGNUM*TJJ ) THEN
                        REC = ONE / XJ
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
                  CALL <a name="SLADIV.438"></a><a href="sladiv.f.html#SLADIV.1">SLADIV</a>( X( J1 ), X( N+J1 ), TMP, Z, SR, SI )
                  X( J1 ) = SR
                  X( N+J1 ) = SI
                  XJ = ABS( X( J1 ) ) + ABS( X( N+J1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale x if necessary to avoid overflow when adding a
</span><span class="comment">*</span><span class="comment">                 multiple of column j1 of T.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( XJ.GT.ONE ) THEN
                     REC = ONE / XJ
                     IF( WORK( J1 ).GT.( BIGNUM-XMAX )*REC ) THEN
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( J1.GT.1 ) THEN
                     CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
                     CALL SAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
     $                           X( N+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span>                     X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 )
                     X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 )
<span class="comment">*</span><span class="comment">
</span>                     XMAX = ZERO
                     DO 50 K = 1, J1 - 1
                        XMAX = MAX( XMAX, ABS( X( K ) )+
     $                         ABS( X( K+N ) ) )
   50                CONTINUE
                  END IF
<span class="comment">*</span><span class="comment">
</span>               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Meet 2 by 2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span>                  D( 1, 1 ) = X( J1 )
                  D( 2, 1 ) = X( J2 )
                  D( 1, 2 ) = X( N+J1 )
                  D( 2, 2 ) = X( N+J2 )
                  CALL <a name="SLALN2.477"></a><a href="slaln2.f.html#SLALN2.1">SLALN2</a>( .FALSE., 2, 2, SMINW, ONE, T( J1, J1 ),
     $                         LDT, ONE, ONE, D, 2, ZERO, -W, V, 2,
     $                         SCALOC, XNORM, IERR )
                  IF( IERR.NE.0 )
     $               INFO = 2
<span class="comment">*</span><span class="comment">
</span>                  IF( SCALOC.NE.ONE ) THEN
                     CALL SSCAL( 2*N, SCALOC, X, 1 )
                     SCALE = SCALOC*SCALE
                  END IF
                  X( J1 ) = V( 1, 1 )
                  X( J2 ) = V( 2, 1 )
                  X( N+J1 ) = V( 1, 2 )
                  X( N+J2 ) = V( 2, 2 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale X(J1), .... to avoid overflow in
</span><span class="comment">*</span><span class="comment">                 updating right hand side.
</span><span class="comment">*</span><span class="comment">
</span>                  XJ = MAX( ABS( V( 1, 1 ) )+ABS( V( 1, 2 ) ),
     $                 ABS( V( 2, 1 ) )+ABS( V( 2, 2 ) ) )
                  IF( XJ.GT.ONE ) THEN
                     REC = ONE / XJ
                     IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
     $                   ( BIGNUM-XMAX )*REC ) THEN
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Update the right-hand side.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( J1.GT.1 ) THEN
                     CALL SAXPY( J1-1, -X( J1 ), T( 1, J1 ), 1, X, 1 )
                     CALL SAXPY( J1-1, -X( J2 ), T( 1, J2 ), 1, X, 1 )
<span class="comment">*</span><span class="comment">
</span>                     CALL SAXPY( J1-1, -X( N+J1 ), T( 1, J1 ), 1,
     $                           X( N+1 ), 1 )
                     CALL SAXPY( J1-1, -X( N+J2 ), T( 1, J2 ), 1,
     $                           X( N+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span>                     X( 1 ) = X( 1 ) + B( J1 )*X( N+J1 ) +
     $                        B( J2 )*X( N+J2 )
                     X( N+1 ) = X( N+1 ) - B( J1 )*X( J1 ) -
     $                          B( J2 )*X( J2 )
<span class="comment">*</span><span class="comment">
</span>                     XMAX = ZERO
                     DO 60 K = 1, J1 - 1
                        XMAX = MAX( ABS( X( K ) )+ABS( X( K+N ) ),
     $                         XMAX )
   60                CONTINUE
                  END IF
<span class="comment">*</span><span class="comment">
</span>               END IF
   70       CONTINUE
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Solve (T + iB)'*(p+iq) = c+id
</span><span class="comment">*</span><span class="comment">
</span>            JNEXT = 1
            DO 80 J = 1, N
               IF( J.LT.JNEXT )
     $            GO TO 80
               J1 = J
               J2 = J
               JNEXT = J + 1
               IF( J.LT.N ) THEN
                  IF( T( J+1, J ).NE.ZERO ) THEN
                     J2 = J + 1
                     JNEXT = J + 2
                  END IF
               END IF
<span class="comment">*</span><span class="comment">
</span>               IF( J1.EQ.J2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 1 by 1 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale if necessary to avoid overflow in forming the
</span><span class="comment">*</span><span class="comment">                 right-hand side element by inner product.
</span><span class="comment">*</span><span class="comment">
</span>                  XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
                  IF( XMAX.GT.ONE ) THEN
                     REC = ONE / XMAX
                     IF( WORK( J1 ).GT.( BIGNUM-XJ )*REC ) THEN
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  X( J1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X, 1 )
                  X( N+J1 ) = X( N+J1 ) - SDOT( J1-1, T( 1, J1 ), 1,
     $                        X( N+1 ), 1 )
                  IF( J1.GT.1 ) THEN
                     X( J1 ) = X( J1 ) - B( J1 )*X( N+1 )
                     X( N+J1 ) = X( N+J1 ) + B( J1 )*X( 1 )
                  END IF
                  XJ = ABS( X( J1 ) ) + ABS( X( J1+N ) )
<span class="comment">*</span><span class="comment">
</span>                  Z = W
                  IF( J1.EQ.1 )
     $               Z = B( 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale if necessary to avoid overflow in
</span><span class="comment">*</span><span class="comment">                 complex division
</span><span class="comment">*</span><span class="comment">
</span>                  TJJ = ABS( T( J1, J1 ) ) + ABS( Z )
                  TMP = T( J1, J1 )
                  IF( TJJ.LT.SMINW ) THEN
                     TMP = SMINW
                     TJJ = SMINW
                     INFO = 1
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  IF( TJJ.LT.ONE ) THEN
                     IF( XJ.GT.BIGNUM*TJJ ) THEN
                        REC = ONE / XJ
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
                  CALL <a name="SLADIV.599"></a><a href="sladiv.f.html#SLADIV.1">SLADIV</a>( X( J1 ), X( N+J1 ), TMP, -Z, SR, SI )
                  X( J1 ) = SR
                  X( J1+N ) = SI
                  XMAX = MAX( ABS( X( J1 ) )+ABS( X( J1+N ) ), XMAX )
<span class="comment">*</span><span class="comment">
</span>               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 2 by 2 diagonal block
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Scale if necessary to avoid overflow in forming the
</span><span class="comment">*</span><span class="comment">                 right-hand side element by inner product.
</span><span class="comment">*</span><span class="comment">
</span>                  XJ = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
     $                 ABS( X( J2 ) )+ABS( X( N+J2 ) ) )
                  IF( XMAX.GT.ONE ) THEN
                     REC = ONE / XMAX
                     IF( MAX( WORK( J1 ), WORK( J2 ) ).GT.
     $                   ( BIGNUM-XJ ) / XMAX ) THEN
                        CALL SSCAL( N2, REC, X, 1 )
                        SCALE = SCALE*REC
                        XMAX = XMAX*REC
                     END IF
                  END IF
<span class="comment">*</span><span class="comment">
</span>                  D( 1, 1 ) = X( J1 ) - SDOT( J1-1, T( 1, J1 ), 1, X,
     $                        1 )
                  D( 2, 1 ) = X( J2 ) - SDOT( J1-1, T( 1, J2 ), 1, X,
     $                        1 )
                  D( 1, 2 ) = X( N+J1 ) - SDOT( J1-1, T( 1, J1 ), 1,
     $                        X( N+1 ), 1 )
                  D( 2, 2 ) = X( N+J2 ) - SDOT( J1-1, T( 1, J2 ), 1,
     $                        X( N+1 ), 1 )
                  D( 1, 1 ) = D( 1, 1 ) - B( J1 )*X( N+1 )
                  D( 2, 1 ) = D( 2, 1 ) - B( J2 )*X( N+1 )
                  D( 1, 2 ) = D( 1, 2 ) + B( J1 )*X( 1 )
                  D( 2, 2 ) = D( 2, 2 ) + B( J2 )*X( 1 )
<span class="comment">*</span><span class="comment">
</span>                  CALL <a name="SLALN2.636"></a><a href="slaln2.f.html#SLALN2.1">SLALN2</a>( .TRUE., 2, 2, SMINW, ONE, T( J1, J1 ),
     $                         LDT, ONE, ONE, D, 2, ZERO, W, V, 2,
     $                         SCALOC, XNORM, IERR )
                  IF( IERR.NE.0 )
     $               INFO = 2
<span class="comment">*</span><span class="comment">
</span>                  IF( SCALOC.NE.ONE ) THEN
                     CALL SSCAL( N2, SCALOC, X, 1 )
                     SCALE = SCALOC*SCALE
                  END IF
                  X( J1 ) = V( 1, 1 )
                  X( J2 ) = V( 2, 1 )
                  X( N+J1 ) = V( 1, 2 )
                  X( N+J2 ) = V( 2, 2 )
                  XMAX = MAX( ABS( X( J1 ) )+ABS( X( N+J1 ) ),
     $                   ABS( X( J2 ) )+ABS( X( N+J2 ) ), XMAX )
<span class="comment">*</span><span class="comment">
</span>               END IF
<span class="comment">*</span><span class="comment">
</span>   80       CONTINUE
<span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SLAQTR.663"></a><a href="slaqtr.f.html#SLAQTR.1">SLAQTR</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
